AD-A257  634 


t 

NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


DTIC 

SELECTE  M 

DECO  1 1992 IJ 

ANALYSIS  OF  CONSTANT  PHASE  CONTOURS  OF 
EVAPORATION  DUCT  MODE  FUNCTIONS 
FOR  WAVEGUIDE  MODE  PROPAGATION 


Jen-Peng  Che 
September  1992 

Thesis  Advisor  Hung-Mou  Lee 


Approved  for  public  release;  distribution  is  unlimited. 


UNCLASSFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


la.  REPORT  SECURITY  CLASSIFICATION 

UNCLASSFIED 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


1b  RESTRICTIVE  MARKINGS 


3.  DISTRIBUTION/AVAILABILITY  OF  REPORT 


2b.  DECLASSIFICATIONAiOWNGRAOING  SCHEDULE 


4.  PERFORMING  ORGANI2ATION  REPORT  NUMBER(S) 


Approved  for  public  release;  distribution  is  unlimited. 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


6a.  NAME  OF  PERFORMING  ORGANIZATION 

Naval  Postgraduate  School 


6c.  ADDRESS  (City.  State,  and  ZIP  Code) 


Monterey,  CA  93943-5000 


6b.  OFFICE  SYMBOL  7a.  NAME  OF  MONITORING  ORGANIZATION 

(If  applicable) 

EC  Naval  Postgraduate  School 


7b.  ADDRESS  (City.  State,  and  ZIP  Code) 


Monterey,  CA  93943-5000 


8a.  NAME  OF  FUNDING/SPONSORING 
ORGANIZATION 


8b.  OFFICE  SYMBOL 
(If  applicable) 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


8c.  ADDRESS  (City.  State,  and  ZIP  Code) 


10.  SOURCE  OF  FUNDING  NUMBERS 


Program  Element  No 


work  unit  Acceuion 
Number 


1 1 .  TITLE  (Include  Security  Qastification) 

ANALYSIS  OF  CONSTANT  PHASE  CONTOURS  OP  EVAPORATION  DUCT  MODE  FUNCTIONS  FOR  WAVEGUIDE  MODE 
PROPAGATION 


12.  PERSONAL  AUTHOR(S)  Jeu-Peug  Che 


13a.  TYPE  OF  REPORT 

Master’s  Thesis 


14.  DATE  OF  REPORT  (year,  month,  day)  1 1 5  PAGE  COUNT 


September  1992 


13b.  TIME  COVERED 
From  To 


16.  SUPPLEMENTARY  NOTATION 

The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official  policy  or  position  of  the  Department  of  Defense  or  the  U.S. 


17.COSATICOOES 

FIELD 

GROUP 

SUBGROUP 

1 8.  SUBJECT  TERMS  (continue  on  reverse  if  necessary  and  identify  by  block  number) 


1 9 .  A8STRACT  (continue  on  reverse  if  necessary  and  identify  by  block  number) 

The  M-Layer  program  tracks  the  constant  phsise  lines  Im{D(q)}  =  0  and  looks  for  their 
intersections  with  the  lines  Re{D(q)}  =  0  for  the  locations  of  me  zeros  of  the  mode  function  D(q). 
These  two  types  of  constant  phase  lines  are  tracked  and  plotted  over  a  search  r^on  which 
contains  modes  having  a  ^ge  attenuation  rate  of  no  more  than  5  dB  per  km.  Several  new 
parameters  for  use  in  mode  search  are  deduced  from  the  results  smd  some  old  ones  are  verified. 
Future  studies  in  waveguide  mode  proparation  theory  pertaining  to  atmospheric  ducts  may 
benefit  from  this  work.  An  improved  mode  sesuxh  strategy  is  also  proposed. 


20^DISTRIBUTION/AVAILABILITY  OF  ABSTRAO 

UNCLASSIFIEOAJNIIMITEO  Q  SAME  AS  REPORT 


22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL 

Hung-Mou  Lee 


DO  FORM  1473, 84  MAR 


21  ABSTRACT  SECURITY  CLASSIFICATION 
UNCLASSFIED 


22b  TELEPHONE  (Include  Area  code) 

(408)646-2846 


83  APR  edition  may  be  used  until  exhausted  SECURITY  CLASSIFICATION 

All  other  editions  are  obsolete  UNCLASSFIED 


22c.  OFFICE  SYMBOL 

EC/Lh 


1 


Approved  for  public  release;  distribution  is  unlimited. 


ANALYSIS  OF  CONSTANT  PHASE  CONTOURS  OF 
EVAPORATION  DUCT  MODE  FUNCTIONS 
FOR  WAVEGUIDE  MODE  PROPAGATION 


by 

Jen-Peng  Che 

Lieutenant  Commander,  Taiwan  Navy 
B.S.,  Chinese  Naval  Academy,  1980 

Submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  ELECTRICAL  ENGINEERING 

from  the 


Author: 


Approved  by: 


NAVAL  POSTGRADUATE  SCHOOL 
September  1992 


Hung-Mou  Lee,  Thesis  Advisor 


Michael  A.  Morgan,  Chairman 
Department  of  Electrical  and  Computer  Engineering 


ABSTRACT 


The  M-Layer  program  tracks  the  constant  phase  lines  Im{D(q)}  =  0  and  looks  for 
their  intersections  with  the  lines  Re{D(q)}  =  0  for  the  locations  of  the  zeros  of  the  mode 
function  D(q).  These  two  types  of  constant  phase  lines  are  tracked  and  plotted  over  a 
search  region  which  contains  modes  having  a  range  attenuation  rate  of  no  more  than  5  dB 
per  km.  Several  new  parameters  for  use  in  mode  search  arc  deduced  from  the  results  and 


some  old  ones  are  verified.  Future  studies  in  waveguide  mode  propagation  theory 
pertaining  to  atmospheric  ducts  may  benefit  from  this  work.  An  improved  mode  search 


strategy  is  also  proposed. 


Accesion  For 

NTIS 

DTIC 

Unann 

Justific 

CRA&I  rV 

TAB  9 

diiticed  □ 

ation 

By 

Distiibution  / 

Availability  Codes 

Dist 

Avail  c 
Spe 

ind/or 

cial 

TABLE  OF  CONTENTS 


L  INTRODUCTION  .  1 

A.  THE  WAVEGUIDE  MODE  THEORY  OF  PROPAGATION .  2 

B.  MODE  SEARCH  PRINCIPLE  .  7 

n.  MODE  FUNCTION  IN  THE  COMPLEX  qj  j/tj  PLANE .  11 

A.  MODE  LOCATIONS  .  11 

B.  PHASE  LINE  TRACKING  .  13 

1.  Downward  Tracking  .  13 

2.  Upward  Tracking  .  18 

in.  ANALYSIS  AND  CONCLUSIONS .  21 

A.  MODE  SEARCH  PARAMETERS  .  21 

B.  MODE  SEARCH  STRATEGY .  24 

1.  Track  Termination .  2< 

2.  Search  Termination .  25 

3.  Track  Duplication  Avoidance .  27 

APPENDIX  A:  MODE  LOCATIONS  IN  THE  COMPLEX  qii/t|  PLANE  .  28 


IV 


APPENDIX  B:  SUBROUTINE  FNDMOD  AND  FZEROX 


49 


APPENDIX  C:  CONSTANT  PHASE  LINE  TRACKING  FLOW  CHARTS .  61 

APPENDIX  D:  CONSTANT  PHASE  LINES  IN  THE  COMPLEX  q^/tj 

PLANE  .  66 

APPENDDC  E:  SEPARATION  OF  Im{D(q))=0  LINES .  75 

APPENDIX  F:  DIFFERENCE  BETWEEN  CONSECUTIVE  ReCq^/tj)  VALUES  84 

LIST  OF  REFERENCES  .  93 

INITIAL  DISTRIBUTION  LIST .  94 


V 


ACKNOWLEDGEMENTS 


I  would  like  to  express  my  sincere  gratitude  to  my  advisor  Dr.  Hung-Mou  Lee  for 
his  professional  guidance,  assistance  and  support  throughout  the  thesis  preparation  period. 
I  thank  my  second  reader  Dr.  Lawrence  J.  Ziomek  for  his  careful  review  and  constructive 
suggestions  to  the  thesis  work. 

Special  thanks  to  my  wife  Mei-Ying  Huang  and  my  son  Yu-Heng  Che  for  their  love 
and  unending  support  which  allowed  me  to  successfully  complete  my  graduate  education. 


VI 


L  INTRODUCTION 


The  M-Layer  program  developed  by  NOSC  (Naval  Ocean  System  Center)  was 
documented  by  Yeoh  [Ref.  1]  at  NFS.  It  was  later  extensively  revised  by  Lee  and  Han 
[Ref.  2]  for  improved  accuracy  and  speed.  The  new  FORTRAN  code  is  now  identified 
as  the  NFS  version  [Ref.  3]  under  the  auspices  of  NRaD  (Naval  Command,  Control  and 
Ocean  Surveillance  Center,  RDT&E  Division),  the  current  name  for  NOSC.  In  this 
version,  the  logical  structures  and  numerical  algorithms  for  following  the  constant  phase 
lines  of  the  mode  function  and  for  computing  the  mode  locations  were  almost  completely 
rewritten.  The  overall  plan  to  first  partition  the  search  region  into  rectangular  areas  and 
then  circulate  around  these  rectangles  to  search  for  the  modes,  i.e.,  the  mode  search 
protocol,  was  left  unchanged.  Lee  and  Han  (Ref.  2]  recommended  that  a  more  direct 
mode  search  protocol  should  be  devised  to  locate  the  modes  more  expediently.  If  possible, 
such  an  approach  should  locate  the  modes  in  the  order  of  ascending  range  attenuation 
rates. 

To  design  a  more  efficient  mode  search  protocol,  a  better  understanding  of  the  mode 
function  beyond  the  known  locations  of  the  modes  is  required.  In  this  thesis,  the 
analytical  property  of  the  mode  function  is  investigated.  Frograms  arc  written  to  track  and 
plot  the  constant  phase  lines  along  which  the  mode  function  is  either  real  or  purely 
imaginary.  From  the  results,  a  new  strategy  to  locate  the  modes  is  recommended. 
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In  the  remainder  of  this  chapter,  the  background  theory  and  some  of  the  notations 
used  in  this  thesis  are  introduced.  The  basic  theory  presented  in  Sections  A  and  B  follows 
Ref.  3  closely.  The  results  of  this  research  are  presented  in  Chapter  11.  In  Chapter  ID,  a 
new  mode  search  protocol  is  proposed  based  on  the  findings  of  this  thesis.  The  relevant 
parameters  for  use  in  this  new  approach  are  also  discussed. 

A.  THE  WAVEGUIDE  MODE  THEORY  OF  PROPAGATION 

Trapping  of  electromagnetic  (EM)  waves  in  the  modes  supported  by  a  duct  is  the 
dominating  factor  in  over-the-horizon  propagation.  The  computer  program  M-Layer 
searches  for  these  modes  and  computes  the  electric  field  from  the  Hertz  vector.  Assume 
a  vertical  electric  dipole 

J  =  4nzfi(x)d(y)8(z-Zj.) 

at  a  height  z=z>p  above  the  surface  of  the  earth.  Following  Freehafer  [Ref.  4],  the  Hertz 
vector  of  the  EM  fields  can  be  written,  in  the  cylindrical  coordinates  (p,  z),  as  a  sum 

of  contributions  from  individual  modes  [Ref.  1]: 

m 

where  is  the  wavenumber  of  the  m-th  mode  and  is  independent  of  the  coordinates; 
Zq-  is  the  height  of  the  transmitter  above  the  ground,  Hq^^^  is  the  Hankel  function  of  the 
second  kind,  which  represents  an  outgoing  wave  in  the  radial  direction  when  the  time 
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dependence  is  adopted,  and  gj^  is  the  height-gain  function  of  the  m-th  mode, 
normalized  so  that  the  integral  of  its  square  over  all  heights  equals  unity  when  either  an 
electric  or  a  magnetic  dipole  source  is  used.  The  height-gain  function  satisfies  the 
equation 


I*  ) 


(2) 


where  k  is  the  free-space  wavenumber,  m(z)  is  the  modified  index  of  refraction,  and 

'J 

m  (z)  is  approximated  with  a  continuous,  piecewise  linear  profile  having  I  layers: 


m\z)=mf+ajiz-z) 


(3) 


In  practice,  m(z)  deviates  only  slightly  from  unity  in  the  troposphere.  The  modified 
refractivity  M(z)  =  [m(z)-l]xlO^,  as  shown  in  Fig.  1,  is  commonly  used.  The  slope  of 
M(z)  is  approximately  o-xlO^/2. 

In  the  i-th  layer,  the  height-gain  function  is  given  by 


where  q^  is  a  dimensionless  linear  function  of  height  z  with  k,  Pj^,  m^,  and  atj  as 
parameters: 


3 


The  functions  and  k2(qjj^)  are  given  in  terms  of  the  Airy  function  Ai: 

and 

A2(0=2(12)*^^i(-g^).  O) 

'J 

For  z  >  Zj^j,  m  (z)  is  extended  continuously  upward  at  a  constant  slope 

—7  —1 

commensurate  with  the  effective  radius  of  the  eanii.  This  slope  equals  2.36x10  (m  ) 
for  the  four-thirds  effective  earth  radius  model.  Thus,  the  height-gain  function  g^^  is  again 
given  by  Eqs.  (4)  through  (7)  for  z  >  Zj^j,  with  the  constant  Cj^j  set  to  to 

represent  an  upward  going  wave  as  z  becomes  large.  Beneath  the  "flattened"  earth  surface 
where  z  <  zj  =  0,  m  (z)  is  deduced  from  the  dielectric  constant  and  the  conductivity  of 
sea  water,  which  ate  set  to  80.8869  and  4.64  (Si/m)  respectively.  Here  the  Hertz  vector 
is  represented  by  a  downward  propagating  plane  wave  in  this  region. 

The  conditions  that  the  tangential  components  of  the  electric  and  magnetic  fields 
are  continuous  across  the  layer  boundaries  determine  the  coefficients  Cj.  Starting  from 
the  top  layer  down,  or  "integration-down"  [Ref.  1],  every  Cj  is  uniquely  determined  by 
at  z  =  Z|^|.  Thus,  a  set  of  all  Cj  coefficients  is  determined  without  considering  the 
boundary  conditions  at  z  =  Zj.  On  the  other  hand,  starting  from  the  lowest  boundary  at 


z  =  Zj  up  to  z  =  Zj,  or  "integration-up",  a  second  set  of  Cj  coefficients  is  determined 
without  applying  the  boundary  conditions  at  z  =  Zj^^.  That  these  two  sets  of  Cj 
coefficients  are  identical  is  a  manifestation  that  is  the  wavenumber  of  a  mode. 
Mathematically,  this  consistency  condition,  sometimes  called  the  guidance  condition 
[Ref.  5],  can  be  expressed  as  the  vanishing  of  a  determinant  D  called  the  mode  function 
[Ref.  6].  The  elements  of  this  determinant  consist  of  linear  combinations  of  kj(qj^p  and 
k2(qmi)  their  derivatives  with  respect  to  height  at  the  boundaries.  The  M-Layer 
program  searches  for  the  zeros  of  the  mode  function  to  find  the  wavenumbers.  Once  a 
wavenumber  is  obtained,  the  height-gain  function  of  the  corresponding  mode  can  be 
computed.  Note  that  the  normalization  condition  on  gj^  and  the  boundary  conditions, 
together  with  the  Cj  coefficients,  determine  the  Bj  coefficients  to  within  an  overall  sign. 
This  sign  need  not  be  resolved  because  only  products  in  the  form  of  gj^(z)gj^(zj)  from 
each  mode  contribute  to  the  Hertz  vector. 

Since  the  mode  function  D  depends  on  the  wavenumber  only  through  the  values 
of  q^^  at  the  layer  boundaries,  it  is  more  convenient  to  consider  D  as  a  function  of  the 
variable  q  given  by 


9  = 


(  k 


\2/ 


v®i; 


«i- 


(8) 


and  to  search  for  the  zeros  of  D(q)  in  the  complex  q  plane.  The  m-th  zero  is  designated 
as  q^^  and  is  called  a  q-eigenvalue  and  Pj^  is  then  deduced  from  q^^  by  inverting  Eq.  (8) 
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for  p.  In  the  NPS  version,  is  ordered  in  ascending  attenuation  rate  of  the  mode,  which 
is  approximately  proportional  to  the  imaginary  part  of 

B.  MODE  SEARCH  PRINCIPLE 

M-Layer  searches  for  modes  which  have  attenuation  rates  below  a  specified  value. 
At  ground  ranges  more  than  several  wavelengths  away,  this  attenuation  rate  is 
proportional  to  the  imaginary  part  of  the  wavenumber  Pj^  as  can  be  seen  from  Eq.  (1). 
In  the  upper  complex  q  plane,  for  values  of  q  such  that 


mj-P/k  is  approximately  proportional  to  q  because  the  absolute  value  of  mj  and,  hence, 
that  of  P/k,  are  both  nearly  unity.  Thus,  the  limit  on  attenuation  rate  imposed  on  the 
imaginary  part  of  P  can  be  translated  approximately  into  an  upper  bound  ^  on  the 
imaginary  part  of  q  which  defines  a  strip  in  the  complex  q  plane  called  the  search  region. 
This  region  is  covered  with  layers  of  identical  square  meshes  whose  sides  are  parallel  to 
the  imaginary  q  axis  and  have  a  length  equivalent  to  an  attenuation  rate  of  1/32  dB  per 

^  Since  the  absorption  is  small  in  air,  the  modified  index  of  refraction  m(z),  and 
hence  a-,  are  considered  as  real  quantities  in  the  following  discussions.  In  the  actual 
FORTRAN  code,  they  are  declared  as  complex  variables. 

2 

This  is  the  default  value  of  the  NPS  version  which  can  be  adjusted  through  editing 
the  variable  "dmesh"  in  an  ASCII  input  file. 
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-1 

kilometer.  The  lower  edges  of  the  meshes  in  the  lowest  layer  sit  along  the  real  q  axis  . 
The  meshes  in  the  layer  just  below  the  top  one  contain  the  upper  bound  of  the  search 
region,  with  the  top  layer  providing  some  allowance  for  numerical  inaccuracy.  The 
program  tests  each  mesh  square  for  the  presence  of  zeros  of  the  mode  function  D(q). 

To  search  through  all  the  meshes,  the  program  first  divides  the  mesh  covering  the 
search  region  into  "contour  rectangles"  with  equally  spaced  vertical  lines  parallel  to  the 
imaginary  q  axis.  These  vertical  lines  contain  the  edges  of  stacked  square  meshes  and  are 
separated  by  a  distance  160  times  the  side  of  a  mesh.  The  search  commences  at  the  top 
left  comer  and  moves  counterclockwise  around  each  "contour  rectangle"  and  begins  with 
the  one  whose  left  edge  is  defined  by 


Rb 


'Jl 


-2/3 


-2xlO-«12e(Mj  -Af^)  «2xlO-V, 


(10) 


where  nij^j^  is  the  minimum  value  of  the  modified  index  of  refraction  profile  and 
is  shown  in  Fig.  1.  The  justification  for  this  choice  to  begin  the  mode  search  is  as 
follows:  From  Eq.  (5)  and  optics,  for  a  wave  to  propagate  freely  in  space,  the  dispersion 
relation  of  the  Maxwell  equations  requires  that  k^m^(z)  >  ^suming  that  both 
quantities  are  real.  Thus,  the  smallest  P  ^  to  support  a  trapped  wave  will  occur  when 


^  The  NOSC  version  extends  the  lower  edge  slightly  below  the  real  q-axis.  This 
causes  problems  in  some  situations  [Ref.  2]. 
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just  exceeds  the  minimum  of  m^(z)^.  It  is  not  surprising  that  an  argument  based 
on  optics  works  within  a  waveguide  mode  formulation  which  is  low  frequency  in 
principle.  Lee  [Ref.  7]  has  demonstrated  that  the  earth-flattening  approximation  actually 
links  Mie’s  low  frequency  oriented  spherical  harmonics  to  Fock’s  high  frequency 
diffraction  theory  through  the  uniform  asymptotics.  The  mere  introduction  of  an  unlimited 
ground  range  into  the  Maxwell  differential  equations  incorporates  the  global  feature  of 
the  radius  of  the  earth  into  the  local  equations. 

After  the  search  over  the  initial  rectangle  is  completed,  the  program  goes  on  to 
search  the  neighboring  rectangle  to  the  left  If  a  specified  number^  of  consecutive 
rectangles  of  decreasing  real  q  values  have  been  searched  without  turning  up  any  zero  of 
D(q),  the  program  changes  direction  and  starts  to  search  the  rectangles  to  the  right  of  the 
initial  "contour  rectangle"  one  by  one,  with  increasing  real  part  of  q.  After  failing 
consecutively  to  locate  any  q-eigenvalue  again  over  the  specified  number  of  rectangles, 
the  program  assumes  that  no  more  zeros  of  D(q)  exists  in  the  search  region.  The  mode 
search  is  considered  complete  and  the  procedure  is  terminated.  If  the  array  for  storing  the 
q-eigenvalues  is  filled  up  before  the  search  is  completed,  the  search  is  terminated  with 
an  error  message. 


^  Note  that  assuming  the  same  real  part,  an  imaginary  component  of  reduces  the 
real  part  of  which  may  enable  the  wave  to  propagate  in  the  z-direction. 

^  This  number  of  consecutive  rectangles  is  an  adjustable  input  variable  "nstop"  in  the 
NFS  version.  The  default  value  is  2. 
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The  search  for  zeros  of  D(q)  makes  use  of  the  fact  that  a  real  valued  function 
changes  sign  when  it  crosses  a  simple  zero.  Since  a  zero  of  the  complex  valued  function 
D(q)  is  where  both  its  real  part  and  imaginary  part  vanish,  a  necessary  condition  for  a 
point  to  be  a  zero  is  that  it  is  the  intersection  of  two  curves  defined  by  Im{D(q) }  =  0 
and  Re{D(q)}  =0.  M- Layer  moves  along  each  side  of  a  "contour  rectangle"  while 
searching  for  a  sign  change  in  Im{D(q)}  across  an  edge  of  a  mesh  bordering  the  side  of 
this  rectangle  to  determine  that  a  line  of  Im{D(q)}  =  0  has  been  encountered.  The  search 
then  follows  this  line  into  the  meshes  within  the  "contour  rectangle",  checking  each  mesh 
to  see  if  a  curve  Re{D(q)}  =  0  enters  the  mesh  being  inspected.  All  these  steps  make  use 
of  the  assumption  that  the  zeros  of  D(q)  are  simple.  Once  both  the  curve  Im{D(q)}  =  0 
and  the  curve  Re{D(q)}  =  0  are  found  to  be  present  within  a  mesh,  the  locations  of  their 
possible  intersections  are  estimated. 

The  rule  for  sign  change  becomes  inconclusive  if  a  zero  of  Im{D(q))  =  0  or  if 
Re{D(q)}  =  0  happens  to  fall  on  a  comer  of  a  mesh,  and  a  remedy  is  required.  In  the 
NFS  version,  whenever  the  real  part  or  the  imaginary  part  of  D(q)  vanishes  on  a  comer 
of  a  mesh,  the  phase  angle  of  D(q)  is  rotated  by  radians.  This  maneuver  effectively 
shifts  q  by  a  small  amount  to  resolve  the  sign  ambiguity. 

To  estimate  the  locations  of  zeros  of  D(q)  within  a  mesh,  D(q)  is  assumed  to  be 
well  approximated  in  this  mesh  by  its  four-term  Taylor  series  expansion.  The  unshifted 
values  of  D(q)  at  the  mesh  comers  determine  this  cubic  polynomial  uniquely.  Cardan’s 
formulas  [Ref.  8]  are  used  to  locate  the  zeros  of  this  polynomial  before  retaining  only 
those  lying  within  the  mesh. 
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n.  MODE  FUNCTION  IN  THE  COMPLEX  q  j  jAj  PLANE 


A.  MODE  LOCATIONS 

M-Layer  searches  for  the  zeros  of  the  mode  function  D(q)  in  the  complex  q  plane 
by  tracking  the  constant  phase  lines  of  D(q)  along  which  the  mode  function  is  either  real 
or  purely  imaginary.  The  zeros  found  are  called  the  q-eigenvalues  as  defined  by  Eq.  (8) 
and  are  denoted  as  These  eigenvalues  are  saved  in  an  ASCII  file  and  are  utilized  later 
for  height-gain  function  computations.  In  order  to  design  a  strategy  for  locating  these 
zeros,  the  mode  locations  are  plotted. 

It  is  clear  from  Eq.  (8)  that  the  variable  q  varies  with  a^,  the  slope  of  m  (z)  in  the 
lowest  (first)  layer.  Since  depends  on  how  the  continuous,  piecewise  linear 
approximation  to  m  (z)  is  made,  while  both  mj  and  are,  in  principle,  dependent  only 
on  the  actual  profile,  q  will  vary  strongly  with  aj  and  fail  to  provide  information  on  the 
wavenumber  of  the  mode  directly.  A  more  suitable  variable  to  use  is,  from  Eq.  (8): 


(11) 


Since  q  is  given  the  variable  name  of  qj  j  and  (k/a  j)^  is  given  the  variable  name  of  tj 
in  the  M-Layer  FORTRAN  code,  this  variable  is  identified  as  qu/tj  throughout  this 
thesis.  Since  m^  is  real  and  both  it  and  P/k  deviate  firom  unity  only  slightly  for  all  cases 
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investigated,  q^/tj  is  approximately  2(mj  -  fi/k).  Thus,  its  imaginaiy  part  is  directly 
proportional  to  -P,  the  attenuation  rate  of  the  mode.  Furthermore,  qj  j/tj  does  not  depend 
explicitly  on  Oj.  Removing  the  factor  from  qj  j  makes  it  possible  to  compare  results 
from  different  evaporation  duct  profiles  meaningfully. 

Plots  of  the  q-eigenvalues  as  individual  points  in  the  complex  plane  are 

included  in  Appendix  A.  For  the  20  meter  duct,  which  will  be  used  as  the  representative 
case  for  discussions  in  this  thesis,  a  line  linking  all  the  eigenvalues  in  the  order  of 
increasing  attenuation  rate  is  shown  in  Fig.  2.  Because  of  the  long  excursion  between 
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many  pairs  of  consecutive  modes  in  the  upper  part  of  this  figiue,  it  appears  that  the 
intuitive  approach  of  starting  with  the  mode  of  the  lowest  attenuation  and  continually 
looking  for  the  mode  with  the  next  higher  attenuation  rate  may  not  be  the  most  efficient 
strategy.  Furthermore,  there  seems  to  be  a  fork  near  the  middle  of  the  figure  which 
suggests  that  there  can  be  more  than  one  constant  phase  line  passing  through  all  the  zeros 
of  the  mode  function.  These  problems  suggest  that,  in  order  to  design  an  efficient  mode 
search  strategy,  the  behavior  of  the  mode  function  in  the  complex  plane  has  to  be 
investigated  more  thoroughly. 

B.  PHASE  LINE  TRACKING 

The  M-Layer  program  follows  a  line  along  which  Im{D(q)}  =  0  until  a  zero  is 
encountered  at  its  intersection  with  another  line  along  which  Re{D(q)}  =  0.  Since  there 
is  no  way  to  determine  the  phase  of  the  constant  phase  line  linking  two  adjacent  zeros  a 
priori,  following  these  two  types  of  somewhat  arbitrary  but  easy  to  compute  phase  lines 
which  have  constant  phases  of  0  or  ti,  and  n/2  or  3n/2,  respectively,  are  still  the  best  way 
to  program  a  computer  to  locate  the  zeros  systematically.  It  is  thus  necessary  to  find  out 
the  actual  distribution  of  these  constant  phase  lines  in  the  complex  q^  jA^  plane. 

1.  Downward  Tracking 

In  the  original  design,  along  the  real  q  axis,  M-Layer  divides  the  search  region 
into  "contour  rectangles,"  each  of  which  spans  160  meshes  horizontally.  From  the 
observed  locations  of  the  modes,  it  is  found  desirable  to  track  and  plot  the  constant  phase 
lines  along  Im{D(q)}  =0  and  Re{D(q)}  =0  over  a  little  more  than  six  "contour 


13 


rectangles",  roughly  three  to  the  right  and  three  to  the  left  of  the  starting  real  q  coordinate 
given  by  Eq.  (10),  which  is  a  variable  called  q^g^^  in  M-Layer,  The  subroutine  FNDMOD 
and  FZEROX  are  modiried  to  track  these  constant  phase  lines.  A  listing  of  these  modified 
routines  for  tracking  the  lines  Im{D(q)}  =  0  beginning  from  the  top  edge  and  moving 
downwards  into  the  search  region  is  included  in  Appendix  B.  Their  flow  charts  are  given 
in  Appendix  C.  Only  minor  changes  are  required  to  track  the  type  of  lines  Re{D(q)}  =  0, 
or  to  track  these  lines  from  the  bottom  of  the  search  region  upwards.  The  program  first 
searches  for  a  sign  change  in  Im{D(q)}  or  in  Re{D(q)}  along  the  top  edge  of  the  search 
region,  starting  at  the  coordinate  -512  mesh  sizes  from  q^g^j,  until  a  distance  of  1024 
mesh  sizes  is  covered.  When  a  sign  change  is  observed,  a  constant  phase  line  is 
recognized  as  passing  through  the  particular  mesh  square  and  the  program  writes  the 
complex  q  value  of  the  lower-left  comer  of  this  mesh  square  into  an  ASCII  file  identified 
by  whether  the  mode  function  is  real  or  imaginary  along  the  line,  and  the  order  this  line 
is  found.  The  program  then  moves  from  the  top  edge  downwards  to  follow  this  constant 
phase  line  until  it  exits  the  search  region.  The  q  values  of  the  lower-left  comer  of  the 
mesh  squares  along  this  phase  line  are  also  written  into  the  same  ASCII  file  to  be  plotted 
later  using  MATLAB/386.  The  constant  phase  lines  for  the  mode  functions  of  the  2,  4, 
6,  8,  10,  20,  30  and  40  meter  evaporation  ducts  are  obtained  and  plotted.  They  are 
included  in  Appendix  D.  Tracking  and  plotting  these  constant  phase  lines  are  extremely 
time  consuming.  The  CPU  time  required  to  track  the  desired  constant  phase  lines  for  each 
duct  using  an  Intel  80486  based  PC  running  at  a  33  MHz  clock  rate  is  listed  in  Table  1. 
The  time  required  to  plot  those  lines  for  each  duct  using  an  Intel  80386  based  PC  running 
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at  a  16  MHz  clock  rate  is  also  listed. 


TABLE  1 

CPU  lune  for  Tracking  and  Plotting  Constant  Phase  Lines 


duct 

tracking 

tracking 

plotting 

height 

Im{D(q)}  =  0 

Re{D(q)}  =  0 

both  lines 

(m) 

(hr:min:sec) 

(hr:min:sec) 

(hr:min:sec) 

2 

1:47:39 

1:46:23 

0:21:28 

4 

2:47:41 

2:48:50 

0:23:17 

6 

3:43:54 

3:45:11 

0:25:42 

8 

4:25:49 

4:23:41 

0:29:34 

10 

6:08:14 

6:11:03 

0:32:25 

20 

8:26:58 

8:24:22 

0:34:21 

30 

9:08:14 

9:06:45 

0:36:49 

40 

8:39:04 

8:37:54 

0:36:53 

Total 

45:07:33 

45:04:09 

4:01:29 

The  constant  phase  line  plot  for  the  mode  function  of  the  20  m  evaporation 
duct  is  given  in  Fig.  3.  The  solid  lines  represent  those  having  Im{D(q)}  =  0;  the  dotted 
lines  represent  those  having  Re{D(q)}  =  0.  Every  intersection  of  these  two  types  of  phase 
lines  is  the  location  of  a  q-eigenvalue  scaled  by  tj,  thus  indicating  the  presence  of  a 
mode.  Note  that  every  solid  line  intersects  with  a  dotted  line  except  for  many  of  those 
running  from  the  top  through  the  bottom  of  the  search  region.  None  of  the  solid  lines 
intersect  with  other  solid  lines,  neither  do  the  dotted  lines.  This  is  seen  clearly  in  Fig.  4, 
which  magnifies  the  lower  center  portion  of  Fig.  3.  One  feature  that  can  be  recognized 
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Constant  phase  lines  in  the  qll/tl  plane  for  20  m  duct 


xl0  5 


immediately  in  all  the  plots  in  Appendix  D  is  that  the  starting  real  q  coordinate  given  by 
Eq.  (10),  which  is  marked  by  an  "x"  along  the  top  edge  of  each  plot,  correlates  very 
closely  with  the  mode  of  ntinimum  attenuation. 

2.  Upward  Tracking 

One  mode  of  low  attenuation  rate  is  missing  from  Fig.  2  compared  to  those 
found  in  Ref.  2.  Several  of  them  are  also  missing  from  those  cases  of  greater  duct  heights 
plotted  in  Appendix  D.  It  is  evident  that  tracking  the  constant  phase  lines  from  the  real 
q^  j/tj  axis  upwards  is  necessary.  For  the  10  meter  and  the  20  meter  ducts.  Fig.  5  shows 
several  constant  phase  lines  starting  out  of  the  lower  limit  of  the  search  region  and 
returning  to  the  lower-half  complex  q^  j/tj  plane.  Fig.  6  shows  the  presence  of  such  lines 
for  the  30  meter  and  the  40  meter  ducts.  These  constant  phase  lines  have  not  been 
observed  in  those  cases  of  lower  duct  heights. 
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Figure  5.  Upward  going  constant  phase  lines  in  the  q  j  |Aj  plane  for  the  10  meter  (top) 
and  20  meter  (bottom)  ducts. 
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xlO-6  Constant  phase  lines  in  the  ql  1/tl  plane  for  30  m  duct 
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figure  6.  Upward  going  constant  phase  lines  in  the  (JjjAj  plane  for  the  30  meter  (top) 
and  40  meter  (bottom)  ducts. 
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m.  ANALYSIS  AND  CONCLUSIONS 


A-  MODE  SEARCH  PARAMETERS 

The  implementation  of  a  mode  search  procedure  requires  several  parameters.  The 
search  region  has  to  be  defined  first.  In  M-Layer,  the  desired  maximum  range  attenuation 
rate  of  the  modes  which  support  wave  propagation  is  treated  as  an  input  parameter  called 
^loss’  conventionally  in  dB  per  km.  This  parameter  determines  the  region  over 
which  the  modes  are  to  be  searched  as  explained  in  Chapter  I.  In  the  complex  q  plane, 
under  the  assumption  that  both  m^  and  p/k  of  interest  are  close  to  unity,  the  search  region 
is  bounded  from  above  by  the  FORTRAN  variable  q^^p  given  by 


(12) 


The  location  at  which  to  start  the  mode  search,  q^ggj.  is  another  parameter 
determined  by  M-Layer.  Based  on  optical  considerations,  the  M-deficiency,  of 

Fig.  1,  which  is  the  amount  of  decrease  of  the  modified  refractivity  from  its  value  on  the 
sea  surface  to  its  minimum  value  in  the  air,  is  a  measure  of  the  capability  of  the  duct  to 
trap  electro-magnetic  waves.  As  explained  also  in  Chapter  I,  this  quantity  determines  the 
location  for  the  program  to  start  searching  for  modes.  The  validity  of  the  optical  argument 
and  the  usefulness  of  Eq.  (10)  have  been  proved  by  this  work,  as  pointed  out  in  Chapter 
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n. 


To  search  for  zeros  of  the  mode  function,  M-Layer  first  covers  the  search  region 
with  identical  mesh  squares.  It  then  follows  the  lines  of  the  type  Im{D(q)}  =  0  through 
each  mesh  square,  checking  for  indications  that  a  curve  Re{D(q)}  =  0  is  entering  the 
same  mesh  so  that  an  intersection  is  possible.  The  term  "mesh  size"  will  denote  the  size 
of  the  edges  of  the  mesh  squares.  It  is  the  size  of  the  step  taken  by  the  program  to 
advance  through  the  search  region.  It  also  determines  the  initial  resolution  of  the  location 
of  a  mode.  The  choice  of  the  mesh  size  should  strike  a  balance  between  the  desire  for  a 
speedy  completion  of  the  search  process  and  the  requirement  in  mode  locating  accuracy. 
As  reported  in  Ref.  2,  for  all  the  evaporation  ducts  considered,  the  choice  of  a  mesh  size 
equivalent  to  an  attenuation  rate  of  1/32  dB  per  km  appears  to  be  optimal,  that  is,  all 
modes  can  be  found  for  all  cases  investigated  in  Ref.  2  without  experiencing 
extraordinarily  long  computation  time.  Consider  this  mesh  size  as  the  default  and  call  it 
q^  in  the  complex  q  plane.  Under  the  same  assumption  for  deducing  Eq.  (12),  the  value 
of  q^/tj  at  9.6  GHz  equals  3.58x10“^.  Figure  7  shows  the  separation  between  two 
adjacent  constant  phase  lines  of  the  type  Im{D(q)}  =  0,  indicated  in  the  figure  as  solid 
lines,  along  the  top  edge  of  the  search  region  for  the  20  meter  duct  The  minimum  of 
these  separations  in  terms  of  q^/tj  is  about  2.8x10  ,  which  corresponds  to  a  spacing 
of  a  little  less  than  eight  mesh  squares  apart  The  minimum  separation  between  adjacent 
Im{D(q)}  =  0  and  Re{D(q)}  =  0  lines  is  thus  about  four  meshes.  Constant  phase  line 
separations  for  all  evaporation  ducts  considered  are  included  in  Appendix  E.  Note  that  the 
minimum  separation  is  almost  constant  for  ducts  higher  than  20  meters.  It  increases 
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slightly  as  duct  height  is  decreased. 

As  the  program  follows  an  Im{D(q)}  =  0  line  into  the  search  region,  a  neighboring 
Re{D(q)}  =  0  line  moves  closer  and  enters  into  the  same  mesh  square.  This  causes  the 
program  to  invoke  the  root  finding  routine  to  determine  the  possible  locations  of  zeros 
of  the  mode  function  within  this  mesh  square.  If  the  mesh  size  is  too  large  and  the 
Re{D(q)}  =  0  line  enters  the  mesh  square  long  before  the  actual  intersection  takes  place, 
the  root  finding  routine  will  be  prematurely  invoked  many  times  and  the  probability  of 
producing  false  modes  is  increased.  This  explains  why  a  larger  mesh  size  sometimes  will 
increase  the  execution  time. 

Given  these  three  parameters  q^Qp,  q^gg^.  and  q^,  a  mode  search  strategy  which  does 
without  the  "contour  rectangles"  is  proposed.  It  should  improve  the  efficiency  of  M-Layer. 

B.  MODE  SEARCH  STRATEGY 

From  the  constant  phase  line  plots  of  Appendix  D,  a  strategy  to  search  for  the  mode 
can  be  drawn:  move  along  the  top  edge  starting  at  Search  first  toward  either  the  left 
or  the  right,  then  reverse  course  to  search  along  the  other  direction.  After  the  search  along 
the  top  edge  is  completed,  search  the  lower  edge  from  one  end  to  the  other.  In  what 
follows,  the  implementation  of  this  strategy  will  be  discussed. 

1.  Track  Termination 

The  tracking  of  an  Im{D(q)}  =  0  line  naturally  ends  if  the  top  edge  or  the 
bottom  edge  of  the  search  region  is  reached.  On  the  other  hand,  the  search  region  is 
unbounded  to  the  right  and  left  of  qjggj-  It  is  convenient  to  retain  the  feature  in  the 
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original  program  to  set  a  limit  on  the  number  of  steps  allowed  to  follow  a  constant  phase 
line.  From  Fig.  3,  this  limit  can  be  set  to  2.5  times  the  number  of  mesh  squares  between 
the  vertical  limits  of  the  search  region.  This  number  equals  2.5x32xajQgg. 

2.  Search  Termination 

M-Layer  has  to  determine  that  no  more  modes  within  the  search  region  is  to 
be  found  and  to  terminate  the  search  for  modes.  When  searching  along  the  top  edge  of 
the  search  region  for  the  constant  phase  lines  Im{D(q))  =  0,  the  separation  between  the 
real  parts  of  the  mode  eigenvalues  is  of  interest  Figure  8  shows  the  separation  in  the  real 
part  of  neighboring  q-eigenvalues  scaled  by  tj,  plotted  against  their  locations  q^/tj  along 
the  real  q^  j/tj  axis  for  the  case  of  the  20  m  duct  Similar  plots  for  all  duct  heights  are 
grouped  in  Appendix  F.  Excluding  the  lowest  point  which  involves  the  mode  obtained  via 
searching  along  the  lower  edge,  the  separation  between  neighboring  modes  is  erratic  with 
an  upward  trend  away  from  the  center  of  the  figure,  which  is  close  to  the  search  starting 
position. 

The  increase  in  distance  between  two  modes  towards  the  end  of  the  range 
searched  makes  it  difficult  to  implement  an  adaptive  mechanism  to  terminate  the  search. 
On  the  other  hand.  Fig.  3  shows  that  almost  all  phase  lines  entering  the  search  region 
from  the  top  that  contain  a  mode  within  this  region  are  bunched  together.  Therefore,  the 
parameter  set  equal  to  four  and  may  be  adjusted,  can  be  used  to  stop  the  search 
along  the  top  edge  when  four  consecutive  constant  phase  lines  tracked  turn  up  no  mode. 

The  search  along  the  real  q  axis  can  be  confined  to  within  the  end  points  of 
the  search  along  the  top  edge. 
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Hgure  8.  OifFerence  between  consecutive  values  (20  m  duct). 
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Re{qnri/tl }  (20  m  duct)  xlO-^ 


3.  Track  Duplication  Avoidance 

When  the  program  searches  step  by  sttp  along  the  top  or  bottom  edges  for 
sign  changes  in  Im{D(q)}  to  start  tracking  the  constant  phase  line,  the  exit  of  a  constant 
phase  line  from  a  previous  track  will  be  encountered.  Entering  into  the  search  region  at 
such  a  location  will  simply  retrace  a  constant  phase  line  which  has  already  been  examined 
for  the  existence  of  a  mode.  Thus,  the  exit  locations  of  the  constant  phase  lines  must  be 
recorded  and  checked  whenever  a  sign  change  in  Im{D(q))  is  found  before  deciding  to 
follow  the  constant  phase  line  into  the  search  region.  Along  the  top  edge,  a  record  of  four 
updated  most  recent  exit  locations  from  the  top  edge  should  be  adequate.  The  locations 
of  exits  from  the  bottom  should  be  kept  for  use  during  the  search  along  the  lower  edge 
of  the  search  region.  For  this  record,  a  dimension  of  256  should  be  adequate  for  the 
current  cases.  But  this  dimension  should  be  adjusted  if  other  types  of  ducts  are  studied. 
A  record  of  four  of  the  most  recent  exit  locations  out  of  the  lower  edge  should  also  be 
kept  and  checked  to  avoid  duplicate  tracking  of  the  constant  phase  lines. 
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APPENDIX  A;  MODE  LOCATIONS  IN  THE  COMPLEX  qjjAj  PLANE 
This  appendix  contains  figures  of  mode  locations  of  evaporation  ducts  of  2,  4,  6, 
8,  10,  20,  30,  and  40  m  heights  plotted  in  the  complex  qj  j/tj  plane  for  easy  comparison. 
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ql  l/tl  (4  meter  duct) 
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;10-6  ql  1/tl  (8  meter  duct) 
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10  ^  qll/tl  (18  meter  duct) 
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ql  1/tl  (24  meter  duct) 


Figure  A.  12  qj  jAj  mode  locations  (24  m  duct). 
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ql  1/tl  (30  meter  duct) 
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ql  1/tl  (32  meter  duct) 


Figure  A.  16  q^/tj  mode  locations  (32  m  duct). 
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ql  1/tl  (38  meter  duct) 


Figure  A.19  qijAj  mode  locations  (38  m  duct). 
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APPENDIX  B:  SUBROUTINE  FNDMOD  AND  FZEROX 


This  appendix  contains  the  listing  of  the  subroutines  FNDMOD  and  FZEROX, 
modified  from  the  NPS  version  of  the  M-Layer  FORTRAN  code  to  track  the  constant 
phase-lines  and  to  locate  the  modes. 
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Line#  Source  Line  Microsoft  FORTRAN  Optimizing  Compiler  Version  5.00 

1  subroutine  fn(lmod(aloss,dmmin.tl,dmesh/ilemj)rmode,qeigen) 

2 

3  c  purpose: 

4  c  This  subroutine  sets  up  the  areas  in  the  complex  ql  1-plane 
Sc  to  search  for  the  zeroes  of  the  modal  function. 

6  c 

7  c  inputs... 

8  c  mxlayr  -  maximum  number  of  layers  allowed  in  refractivity 

9  c  profile 

10  c  nzlayr  -  actual  number  of  layers  in  refractivity  profile 

11  c  aloss  -  maximum  attenuation  rate  (in  db/km)  of  modes 

12  c  to  be  found 

13  c  dmmin  -  minimum  of  zim(j)-zim(l) 

14  c  tl  -  koal23 

15  c  dmesh  -  initial  mesh  size  divisor 

16  c 

17  c  outputs... 

18  c  qeigen  -  complex  array  containing  the  zeros  of  the  modal 

19  c  function 

20  c  nrmode  -  actual  number  of  modes  found 

21  c 

22  c  subroutines  called... 

23  c  fzerox 

24  c  findfx 

25  c  common  block  areas... 

26  c  coml 

27  g********** 

28  c 

29  implicit  double  precision  (a-h.o-z) 

30  complex*  16  ctemp,t  1  ,qeigen.zeros,fx  l,fx2 

3 1  parameter(c201og^.68588963806504d0,one=(  1  .d0,0.d0), 

32  $  Step0=l024.d0.step0h=step0/2.d0) 

33  character*3  filemjdine 

34  character*40  rfile 

35  c 

37  c  qeigen  •  complex  array  containing  all  the  zeros  of  the 

38  c  modal  function  found 

39  c  zoos  -  complex  array  containing  the  zeros  of  the  modal 

40  c  function  found  in  the  current  rectangular  region 

41  c  of  the  complex  ql  1-plane 

42  c 

44  c  use  include  file  for  parameters  of 

45  c  mxlayr  max  #  layers 

46  c  mxmode  max  #  modes 

47  c 

48  Sinclude:  ’mlaparm.inc’ 

*****  Begin  listing  of:  mlaparm.inc 
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1  c 

2  c  include  file  to  deHne  the 

3  c  maximum  #  of  layers  (mxlayr) 

4  c  maximum  #  of  modes  (mxmode) 

5  c 

6  parameter  (mxlayrsSS  ) 

7  parameter  (mxmode=127 ) 

*****  End  listing  oft  mlaparm.inc 

49 

50  dimension  qeigen(mxmode),zeios(2*mxmode-«-l) 

51  c***** 

52  common  /coml/waveno 

53  c 

^  ^^im********************mm***m**********:****m************m****m********* 

55  c  rl  -  value  of  mode  search  on  the  real  part  of  qll  at  the 

56  c  left  edge  in  tmesh  units. 

57  c  t2  -  value  of  mode  search  on  the  real  part  of  ql  1  at  the 

58  c  right  edge  in  tmesh  units. 

59  c  hot  -  value  of  the  imaginary  part  of  ql  1  at  the  bottom 

60  c  edge  (this  is  set  to  0). 

61  c  tO  -  value  of  the  imaginary  part  of  ql  1  at  the  top  edge 

62  c  in  tmesh  units. 

63  c  step  •  size  of  search  areas. 

65  c 

66  c 

68  c  start  of  executable  statements 

71c  set  up  search  areas  for  finding  modes  and  solve  for  modes, 

72  c  calculate  approximate  value  for  re(q  11)  where  modes  start 

73  c  occurring 

^iti*mm*****************************m**m***m*m***********m**************mm 

75  C 

76  nrmode=0 

77  c 

78  reconssdreal(tl) 

79  rstarts-2.0d-6*dmmin 

80  qtest=rstart*recons 

81 

82  ttop  1  =(2.0d-3/(waveno*c201og))*recons 

83  ttopsaloss^ttopl 

84  tmesh=ttopl/dmesh 

85  if(tmesh  .gt.  0.1  dO)  tmesh=1.0d-l 

86  if(tmesh  .gt.  l.Od-3)  then 

87  tol=1.0d-4 

88  else 

89  tol=tmesh*0.1d0 

90  end  if 
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91  c 

92  c***** 

93  write(*,1002) 

94  write(16.1(X)2) 

95  c 

96  write(*,1003)  tmesh.tol 

97  write(16,1003)  tmesh,tol 

98  c 

99  tO=dnim(tt(v/lmesh)+l.dO 

100  bot=0.d0 

101  c***** 

102  c 

103  rl=dnint(qtest/tmesh)-stepOh 

104  rl=rl 

105  rr=rl+stepO 

106  caU  findfx(rl,tOJxUl,yl.tmesh) 

107  il=inl(rl) 

108  in=int(step0)+int(rl) 

109  nlinesO 

110  50  continue 

111  do  100  nn=il.in 

112  r2=rl+l.d0 

113  caU  findfx(t2.t0.fx2,x2,y2,tinesh) 

114  if  (yl*y2  .It.  0.d0)  then 

115  r=fl*tmesh 

116  write  (*.5555)  r 

1 17  c  write  (16,5555)  r 

118  5555  foimat  (/5x.’i»  \d28.16) 

1 19  go  to  400 

120  endif 

121  rl*r2 

122  fxl=fx2 

123  xl*x2 

124  yi=y2 

125  100  continue 

126  return 

127  c 

128  400  continue 

129  nrold^nimode 

130  nlinesnline+l 

131  ncl(X)sint(nline/100) 

132  ncl0=int((nline-ncl00*100)/10) 

133  ncl=nline-ncl00*l(X)-ncl0*10 

1 34  kline=char(48+nc  1 00)//char(48+nc  1 0)//char(48+nc  1 ) 

135  rfile=’e:\ans5Sr7/filetn//kline//*.dat’ 

136  c****** 

137  open  (unit=25,fiIe=rfile,statuss’unknown’) 

138 

139  call  fzerox(t0 jl jl IJx  1  ,x  1  .y  1  /x2,x2.y2,zeros.nroldjimew, 

140  $  tmesh) 
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141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 


close  (25) 
c 

nnnode=ninew 

newO=nnno(ie-nrold 

c 

c  mode  search  completed 
c  order  zeros  found  by  order  of  increasing  real  part. 

^*««*«*»*4<***«4i4i«*4i*4i******»««««**«*««*«««4.«««»********«««*««**4>«**»* 

C 

if(nrmode  .gt  1)  then 
jkend=nnnode- 1 

do  420  jk=l  Jkend 
nend=nrmode-jk 


do  410  ja=l.nend 
nij=nrmode-ja 
nijl=nij+l 

if(dimag(zeros(nrj  1  )).lL 
$  dimag(zeros(nij)))  then 
ctemp=zeros(nijl) 
zerosCnij  1  )=zeros(nij) 
zeros(nij)sctemp 
end  if 

410  continue 
420  continue 
end  if 


^«««*«<»»««*«*««**4i*»*««**4i*****««***<|i*4i«****»****«****««««*4>*****«**«* 


c  the  possibility  exists  that  duplicate  (within  the  tolerance  'toP) 
c  zeros  of  the  modal  equation  will  be  found,  eliminate  these 
c  duplicate  zeros. 

C 


jkflagsO 

jkendsnrmode-1 


do  240  jksl  jkend 
jklsjk+1 
ctempszeros(jkl) 

chksqscdabs(zeros(jk-jkflag)-ctemp) 
if(chksq  .It.  tol)  then 
jkflagsjkf1ag+l 
go  to  240 
end  if 

zeros(jk  1  -jkflag)=ctemp 
240  continue 

nimodesnrmode-jkflag 
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191  c 

192  c******************************************* 

193  c  Store  the  zeros  as  the  eigenvalues. 

194  ^**************************************’^*** 

195  c 

196  nnnodesminO(nnnode.mxmode) 

197 

198  do  430  jksljinnode 

199  qeigen(jk)szeros(jk) 

200  430  continue 

201  c 

202  il=:int(rll)+l 

203  if  (il  .It.  in)  then 

204  rl=rll+l.d0 

205  call  find£x(rl,t0,fxl.xl,yl.tmesh) 

206  go  to  50 

207  endif 

208  return 

209  c 

211  c  format  statements 


213  1000  format(/5x,’searching  for  zeroes  in  this  areas  are*. 

214  $  ’  defined  by:75x,’istart=  ’4lQ/5x,’itop=  ‘JIO, 

215  $  5x,’ ibot* ’411/5x,’ilefla  MlO^x/iright*  \il0) 

216 

217  1001  format(5xi4,’  new  zeroes  found  in  this  area.'/) 

218 

219  1002  fonnal(//5x.’*******  start  search  for  modal  eigenvalues’, 

220  $  '  ••****••) 

221 

222  1003  fcKmat(/5x.’tmesh*  ’.dl5.5,5x,’tol=  ’,dl5.5/) 

223 

224  end 

225  c 

226  c 

228  £*•****•*••**••*****•**••*•*****•••***♦•***•*•••*•*••*••• 

229  c 

230  subroutine  fzerox(t0jlj’ll,fxl,xl,yl4x2,x2.y2,zerosjiioid, 

231  $  nmew,tmsh0) 

232  c 

233  c***** 

234  c  fzerox  is  a  routine  for  finding  the  zeroes  of  a  complex  function,  f, 

235  c  which  lie  within  a  specified  rectangular  region  of  the 

236  c  complex  qll  plane,  assuming  that  the  function  has  only 

237  c  simple  zeroes  over  this  rectangle. 

238  c 

239  c  parameters  specifying  the  search  rectangle: 

240  c  tmesh  •  set  equal  to  about  half  the  average  spacing  between 
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241  c  zeroes  within  the  rectangle.  A  smaller  value  may  be  used 

242  c  as  a  safety  measure,  but  too  small  a  value  will  result 

243  c  in  excessively  long  run  time. 

244  c  zeros  •  output  list  of  (complex)  values  of  ql  1  at  which 

245  c  zeroes  are  found. 

246  c  nmew  juold'  the  number  of  zeroes  found 

247  c 

248  c  subroutines  called- 

249  c  findfx 

250  c  roots 

251  c***** 

252  implicit  double  precision  (a-h,o-z) 

253  complex*16  fl0,f01,fll,fxl,fx2,fx00Jxl0Jx01.fxll, 

254  +  one.sol.zeros 

255  parameter(one=(  1  .d0,0.d0)) 

256  Sinclude:  ’mhq)aim.inc’ 

*****  Begin  listing  of;  mlapann.inc 

1  c 

2  c  include  file  to  define  the 

3  c  maximum  #  of  layers  (mxlayr) 

4  c  maximum  #  of  modes  (mxmode) 

5  c 

6  parameter  (mxlayr=35  ) 

7  parameter  (mxmode=127 ) 

*****  End  listing  of:  mlaparm.inc 

257  dimension  sol(3),theta(2).zeros(2*mxmode+l)jline(2,1024) 

258  c 

259  commonAmccom/tmesh 

260  c***** 

261  c 

262  npointsO 

263  nmew=nrold 

264  tmesh=tmshO 

265  bot=0.d0 

266  tst0-l.d0 

267  rll=rl 

268  rsrll 

269  c***** 

270  c 

271  fx01=fxl 

272  fxOliaxl 

273  fxOlisyl 

274  fxll=fx2 

275  fxllrsx2 

276  fxlli=y2 

277  go  to  82 

278  c 

279  c***** 

280  c  enter  mesh  square  from  left  side  or  exit  rectangle  at  right  edge. 

281  c 


55 


282  20  i^r+LdO 

283  21  fx01»fxll 

284  fxOlPsfxllr 

285  fx01i=fxUi 

286  fxOOsfxlO 

287  fxOOi^fxlOr 

288  fxOOi»fxlOi 

289  22  continue 

290  caU  find^(r+l.d04-»-l.d0Jxl Uxl  Irjfxl  li,tmesh) 

291  call  nndfx(rf  I.d0,tjxl0ixl0r,fxl0i.tmesh) 

292  c******* 

293  c  Determine  the  edge  of  exit  of  un(0=0  from  current  mesh. 

294  edgeit=fx01i'*fxlli 

295  edgeib=fx00i*fxl0i 

296  if  (edgeib  .gt.  O.dO)  then 

297  c  lm(0»0  goes  through  the  01  to  10  line. 

298  if  (edgeit  .gt.  0.d0)  then 

299  c  Im(0s0  goes  through  the  10  to  11  edge  (edge  1). 

300  lout=l 

301  else 

302  c  Im(f)sO  goes  through  the  01  to  11  edge  (edge  2) 

303  louts2 

304  end  if 

305  else 

306  c  lm(0=^  goes  through  the  00  to  10  edge  (edge  4) 

307  louts4 

308  if  (edgeit  .IL  0.d0)  then 

309  c  lm(f)«0  also  nins  through  01  to  11  and  10  to  11  edges. 

310  c  Store  crossing  location  and  in/out  information. 

311  knot34sknot34+l 

312  c  loc34r(knot34)=Ir 

313  c  Ioc34i(knot34)sli 

314  end  if 

315  end  if 

316  c******* 

317  go  to  85 

318  c***** 

319  c  enter  mesh  square  from  bottom  side  or  exit  rectangle  at  top  edge. 

320  40  t=t+l.dO 

321  41  fx00=fx01 

322  fxOOi^fxOlr 

323  fxOOisfxOli 

324  fxlO=fxll 

325  fxlOt*=fxllr 

326  fxlOisfxlli 

327  42  continue 

328  call  findfxfr.t'f  l.dOi'xOli'xOlr.fxOli.tmesh) 

329  call  rmdfx(r-t-l.d0,t-«-l.d0.fxll.fxllr4^xlli,tmesh) 

330  c******* 

331  c  Determine  the  edge  of  exit  of  im(0«0  from  current  mesh. 
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332  edgeil=fx00i*fx01i 

333  edgeirsfxlOi*fxlli 

334  if  (edgeir  .gt.  O.dO)  then 

335  c  lm(0=0  goes  through  the  00  to  1 1  line. 

336  if  (edgeil  .gt.  0.d0)  then 

337  c  lm(f>=0  goes  through  the  01  to  11  edge  (edge  2) 

338  louts2 

339  else 

340  c  lm(0=0  goes  through  the  (X)  to  01  edge  (edge  3). 

341  lout=3 

342  end  if 

343  else 

344  c  lm(0=0  goes  through  the  10  to  1 1  edge  (edge  1) 

345  louts  1 

346  if  (edgeil  .IL  0.d0)  then 

347  c  lm(0=0  also  runs  through  00  to  01  and  01  to  11  edges. 

348  c  Store  crossing  location  and  in/out  information. 

349  knot41=knot41+l 

350  c  loc41r(knot41)=lr 

351  c  loc41i(knot41)=ti 

352  end  if 

353  end  if 

354  €**••*•• 

355  go  to  85 

356  c***** 

357  c  enter  mesh  square  from  right  side  or  exit  rectangle  at  left  edge. 

358 

359  60  r»r-I.d0 

360  61  fxllsfxOl 

361  fxlli^fxOlr 

362  fxllisfxOli 

363  fxlOsfxOO 

364  fxlOrsfxOOr 

365  fxlOisfxOOi 

366  65  continue 

367  call  findfx(r.t+l.d0ix01ix01r,fx01i,tmesh) 

368  call  nndfx(r,t,fx(X).fxOOrix(X)i,tmesh) 

369  c******* 

370  c  Determine  the  edge  of  exit  of  im(0=0  from  current  mesh. 

371  edgeit=fx01i*fxlli 

372  edgeib=fxOOi*fxlOi 

373  if  (edgeit  .gt.  0.d0)  then 

374  c  lm(0=0  goes  through  the  01  to  10  line. 

375  if  (edgeib  .gt.  O.dO)  then 

376  c  Im(f)sO  goes  through  the  00  to  01  edge  (edge  3). 

377  lout=3 

378  else 

379  c  lm(0=0  goes  through  the  00  to  10  edge  (edge  4) 

380  lout=4 

381  end  if 
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382  else 

383  c  lin(f)>=0  goes  through  the  01  to  11  edge  (edge  2) 

384  lout=2 

385  if  (edge*  .It.  0.d0)  then 

386  c  Iin(f)s0  also  runs  through  (X)  to  10  and  00  to  01  edges. 

387  c  Store  crossing  location  and  in/out  information. 

388  knotl2=knotl2+l 

389  c  locl2r(knotl2)=Ir 

390  c  locl2i(knotl2)=li 

391  end  if 

392  end  if 

393  c******* 

394  go  to  85 

395  c***** 

396  c  enter  mesh  square  top  side  or  exit  rectangle  at  bottom  edge. 

397  c 

398  80  t=t-l.d0 

399  81  fx01=fx00 

400  fxOlrsfxOOr 

401  fxOlisfxOOi 

402  fxll=fel0 

403  fxllr=fxl0r 

404  fxlli:sfxl0i 

405  82  continue 

406  call  findfx(r.t,fx00,fx00rJx00i,tmesh) 

407  call  findfx(r+l.d0.tixl0ixl0r,fxl0i,tmesh) 

408  c 

409  c******* 

410  c  Determine  the  edge  of  exit  of  im(f)=0  from  current  mesh. 

411  c 

412  83  continue 

413  edgeil=fx00i*fx01i 

414  edgeir=fxl0i*fxlli 

415  if  (edged  .gt.  O.dO)  then 

416  c  lm(f)=0  goes  through  the  00  to  11  line. 

417  if  (edgeir  .gt.  O.dO)  then 

418  c  Im(f)s0  goes  through  the  00  to  10  edge  (edge  4) 

419  lout»4 

420  else 

421  c  lm(f)=0  goes  through  the  10  to  11  edge  (edge  1). 

422  lout*l 

423  end  if 

424  else 

425  c  lm(f)=0  goes  through  the  00  to  01  edge  (edge  3) 

426  lout*3 

427  if  (edgeir  .It.  0.d0)  then 

428  c  Im(f)s0  also  runs  through  00  to  10  and  10  to  1 1  edges. 

429  c  Store  crossing  location  and  in/out  information. 

430  knot23*knot23+l 

431  c  loc23r(knot23)=lr 
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432  c  Ioc23i(knot23)=li 

433  end  if 

434  end  if 

436  c  Test  for  there  being  at  least  one  ie(0=0  line  entering  and 

437  c  leaving  the  mesh  square. 

439  c 

440  8S  continue 

441  npointsnpoint-fl 

442  rn=r*tmesh 

443  tth=t*tmesh 

444  rline(ljipoint)=rrr 

445  rline(2jipoint)=ttt 

446  if  (npoint  .eq.  1024)  lout=5 

447  c 

448  if((t  .IL  bot)  .or.  (t  .gt.  tO))  go  to  100 

449 

450  if  ((fx00r*fxl0r  .gL  O.dO)  .and.  (fx01r*fxllr  .gt.  0.d0) 

451  +  .and.  (fx00r*fx01r  .gt  O.dO))  go  to  (20,40,60.80.100) 

452  +  lout 

453  c******'*"''*'*'**’*"*''"*"'"''''"'"***'*"*'*"''****'*'*’*'**********************''"'"*'**'''****’*"*”*'''*** 

454  c  Computate  the  values  of  the  modal  function  at  the  comers  of 

455  c  a  mesh  square  to  determine  its  Taylor  series  to  the  3rd  order 

456  c  for  estimating  its  root  locations. 

458  c 

459  c  f00=one 

460  fl0=cdexp(fxl0-fx00)-one 

461  f01=cdexp(fx01-fx00)-one 

462  f  1  l=cdexp<fx  1 1  -fx00)-one 

463  c 

464  Q***********  estimate  locations  of  zeroes  by  radicals***************** 

465  c 

466  call  roots(f lOJO  1/11  .soloirsol) 

467  c 

468  do  90  n=l,nrsol 

469  ureal  =  dreal(sol(n)) 

470  uimag  =  dimag(sol(n)) 

471  if  (ureal  .It.  O.dO  .or.  ureal  .gt.  l.OdO)  go  to  90 

472  if  (uimag  .It.  0.d0  .or.  uimag  .gt.  l.OdO)  go  to  90 

473  92  theta(l)=(r-t-ureal)*tmesh 

474  theta(2)=(t+uimag)*Unesh 

475  nmew  *  nmew+1 

476  zeros(nmew)3:dcmpbi(theta(  1  ),theta(2)) 

477  90  continue 

478  c 

479  c******************************************* 

480  c  continue  following  the  phase  line 

481  c*****’^*'*'****’^****'*'***’^'*'************’*'******* 
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482  c 

483  95  continue 

484  go  to  (20.40.60.80.100)  lout 

485  c** 

486  100  continue 

487  c****** 

488  write  (25.8888)  ((rline(ij).  i=l.  2).  j=l.  npoint) 

489  8888  format  (2el5.5) 

490  c 

491  return 

492  end 
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APPENDIX  C:  CONSTANT  PHASE  LINE  TRACKING  FLOW  CHARTS 


nrmode  -  0 
recons  -  dreal(t1) 
qtest“-2x10^(-6)  x  dmmin  x  recons 
ttop1-((2x10"(-3)  /  kx  20log(e))  x 
recons 

top-  alossxttopi 
tmesh  -  ttop1  /  dmesh 


100 


DO  100 
nn=l1  ,in 
r2*r1+1.d0 


CONTINUE 


CALL  FINDFX 


r1  -r2 
fx1  -fx2 
x1  -x2 
y1-y2 


N/y1  xy2 

“<L  <  0 


r-r1  X  tmesh 


RETURN 


WRITE  ( r ) 


4001 
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nrold  -  nrmode 
niine  -  nline  + 1 
ndOO  -  int(niine/100) 
nc10  -int((nllne-nc10(n00)/10) 
nc1  -  nline  -  nc100*100  -  nc10*10 
kline»char(48+nc1 00)//char(48+nc1 0 
)//char(48+nc1) 

rfile-'e:\ml.che\ans\r‘//filem//klln0//.daf 


Open(unit“25,  file-rfilo,  status* 
'unknown') 


CALL  FZEROX 


CLOSE ( 25 ) 


nrmode  -nmew 
newO  -  nrmode  -  nrold 
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APPENDIX  D:  CONSTANT  PHASE  LINES  IN  THE  COMPLEX  q^Aj  PLANE 
This  appendix  contains  plots  of  the  constant  phase  lines  Im{D(q))  =0  and 
Re{D(q)}  =  0  initiating  from  the  top  edge  of  the  search  region  in  the  complex  q^  jAj 
plane  for  the  evaporation  ducts  of  2,  4,  6,  8,  10,  20,  30,  and  40  m  heights. 
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Constant  phase  lines  in  the  ql  I/tl  plane  fur  2  in  duct 
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Xl0  5 


Constant  phase  lines  in  the  ql  1/tl  plane  for  4  m  duct 


xIO-5 


Constant  phase  lines  in  the  ql  I/tl  plane  for  8  m  duct 


Figure  D.4  Constant  phase  lines  in  the  q|  jAi  plane  (8  m  duct). 
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xl0  5 


Constant  phase  lines  in  the  ql  I/tl  plane  for  10  in  duct 


xl0  5 


Constant  phase  lines  in  the  ql  1/tl  plane  for  20  m  duct 
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xIO-5 


APPENDIX  E:  SEPARATION  OF  Im{D(q)}=0  LINES 


This  appendix  contains  figures  of  the  separation  of  Ini{D(q)}=0  lines  along  the  top 
edge  of  the  search  region.  The  2,  4,  6,  8,  10,  20,  30  and  40  meter  evaporation  ducts  are 
included. 
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9-0 1 


(n/nb)3H  JO  33U3j3ijia 

Figure  E.1  Separation  of  Iin{D(q))=0  lines  along  the  top  edge  of  the  search  region 
(2  m  duct). 


76 


Phase-Line  Number  (2  m  duct) 


xlO-6 


(lVllb)3H  JO  3DU3J9JJia 

Figure  E.2  Separation  of  Im(D(q))^  lines  along  the  top  edge  of  the  search  region 
(4  m  duct). 
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Phase-Line  Number  (4  m  duct) 


xlO-7 


Figure  E3  Sq}aradon  of  Iin{D(q))=0  lines  along  the  top  edge  of  the  search  region 
(6  m  duct). 
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Phase-Line  Number  (6  m  duct) 


xlO-7 


Figure  E4  Separation  of  Im(EKq)}=0  lines  along  the  top  edge  of  the  search  region 


(8  m  duct). 
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Phase-Line  Number  (8  m  duct) 


xlO-7 


Figure  E.5  Separation  of  Ini{D(q))=0  lines  along  the  top  edge  of  the  search  region 
(10  m  duct). 
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Phase-Line  Number  (10  m  duct) 


Xl0  7 


(lVllb)3^iJ0  30U3J3Jj!a 

Figure  E.6  Sq>aration  of  I]ii{D(q))=0  lines  along  the  top  edge  of  the  search  region 
(20  m  duct). 
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xlO-7 


Rgure  E.7  Separation  of  Ini{D(<l))=^  along  the  top  edge  of  the  search  region 
(30  m  duct). 
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Phase-Line  Number  (30  m  duct) 


xlO-7 


Figure  E.8  Separation  of  Ini{D(q))=0  lines  along  the  top  edge  of  the  search  region 
(40  m  duct). 
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Phase-Line  Number  (40  m  duel) 


APPENDIX  F:  DIFFERENCE  BETWEEN  CONSECUTIVE  ReCq^j/tj)  VALUES 
Separation  in  real  part  of  neighboring  q-eigenvalues  scaled  by  tj  is  plotted  in  this 
appendix  against  the  eigenvalue  locations  along  the  real  qjjAj  axis. 
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XlO-6 


Figure  F.l  Difference  between  consecutive  Re(qjj^j)  values  (2  m  duct). 
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Re  { qm/t  1 }  (2  m  duct)  x  1 0'^ 


xlO-6 


Figure  F.2  Difference  between  consecutive  values  (4  m  dua). 
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Re{qm/ll )  (4  m  duct)  xlO-5 


Xl0  6 


OOVO'^CS  —  OO'O'^CSO' 

«  ^  ^  dodo 


sanjB^  3Apno3suo3  uasMisg  wusjsjjiQ 

Figure  F.3  Difference  between  consecutive  Rc(qj^^j)  values  (6  m  duct). 
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Re  { qm/l  1 }  (6  m  duct)  x  1 0-^ 


xlO-7 


S3npjy\  (lJ/Ulb)3>J  3AIjn09SU03  U93MJ3a  3DU3J3jJIQ 
Figure  F.4  Difference  between  consecutive  Re(q^])  values  (8  m  duct). 
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Re{qm/tl )  (8  m  duct)  xlO' 


xlO-6 


Figure  F,5  Difference  between  consecutive  Re(qj^^j)  values  (10  m  duct). 


89 


Re(qm/ll )  (10  m  duct)  xlO^ 


xlO-7 
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Figure  F.6  Difference  between  consecutive  ReCq^^j)  values  (20  m  duct). 
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Re|qm/ll )  (20  m  duct)  xlO'5 


xlO-7 
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s3niBy\  3Aiino3Suo3  souajajjiQ 

Figure  F.7  Difference  between  consecutive  ReCq^jj/tj)  values  (30  m  duct). 
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Re{qm/ll )  (30  m  duct)  xlO-^ 
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(n/iub)3^  3AiinD3suo3  uMMpg  aouajsjjiQ 


Figure  F.8  Difference  between  consecutive  Refqj^j)  values  (40  m  duct). 
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